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SUMMARY 


A Monte Carlo procedure has been developed to simulate turbu- 
lent boundary layer wall pressure fluctuations. The approach 
utilizes much of the newly available conditional sampling infor- 
mation to construct the required distribution functions. Various 
disturbance wave forms have been examined, as well as the effect 
of frequency-dependent decay. Good agreement between the simulation 
and experimental data has been achieved for root mean square 
pressure level, power spectrum, and space-time correlation. 

SYMBOLS 


f frequency (Hz) 

n random number 

p pressure 

P rmg root mean square pressure 

R _ space-time correlation 

PP 

S, Strouhal number (oj 6*/U ) 

t 00 

T integration time interval 

t time 

t reference time 

o 

free stream velocity 



u c convection velocity 

u friction velocity 

T 

x coordinate in direction of flow 

Xp development length 

y transverse coordinate 

y frequency adjustment parameter 

6* displacement thickness 

v kinematic viscosity 

t dimensionless time (U^t/S* ) 

x wall shear stress 

w 

a) radian frequency 


I . INTRODUCTION 

The purpose of this report is to present a procedure for simu- 
lating wall pressure fluctuations beneath a low speed, fully turbu- 
lent boundary layer. A simulation was desired because it could 
supply local pressure fluctuation data with arbitrary spatial and 
temporal resolution for a variety of velocity and boundary layer 
thickness conditions . The arbitrary resolution feature permits 
direct coupling between simulated boundary layer pressure fluctua- 
tions and numerical structural analysis codes, resulting in straight- 
forward prediction of local surface motion produced by the simulated 
turbulent boundary layer. This work represents part of a compliant 
wall drag reduction study undertaken at NASA/Langley Research 
Center. 

Bushnell, Hefner, and Ash (ref. 1) have discussed a possible 
mechanism for reducing skin friction drag by means of a compliant 
surface. Based on that discussion, it is apparent that such an 
effect can exist only if the surface motion is compatible with 
the turbulent flow structure. They have suggested that the 
violent turbulent "bursts” near the wall, as discussed by Of fen 
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and Kline (ref. 2), can be inhibited by a suitable periodic surface 
motion. In addition to the connection between surface motion and 
inner turbulent flow structure, it is known that the wall pressure 
fluctuations are produced primarily by large vorticular structures 
in the outer portion of the boundary layer (see Praturi, ref. 3, 
for a description of the outer flow structures). Consequently, 
the surface motion is driven by disturbances emanating from the 
outer part of the boundary layer, and is required to interact 
with flow structures in the inner region. These conditions put 
stringent requirements on wall pressure simulation and surface 
motion calculations. Ash and Balasubramanian (ref. 4) have dis- 
cussed some of those problems. It should be evident that the 
simulation of turbulent wall pressure fluctuations must be connected 
as closely as possible to the physics of the flow. Furthermore, 
the simulated turbulent structure must result in pressure fluctua- 
tion statistics that reproduce accepted experimental data. 

Willmarth (ref. 5) has summarized most of the experimental 
measurements of turbulent wall pressure statistics, and only those 
measurements which were needed for the simulation are mentioned 
here. Willmarth and Woolridge (ref. 6) and Bull (ref. 7) have 
presented detailed measurements of wall pressure energy spectra 
and space-time correlations. Both investigations found that the 
pressure energy spectrum peaked at a Strouhal number (S, = u> 6*/U ) 

L 00 

of about 0.2, that the root mean square pressure was related to 
the wall shear stress (P rms ^2.5 to 2.8 x w ), and that pressure 
fluctuations were convected downstream with varying speeds. They 
were uncertain whether the pressure disturbances with different 
sizes (or frequencies) were convected with different velocities, 
or whether disturbances actually changed speeds, or both. Both 
investigations did reveal more rapid decay in the correlation 
of high frequency pressure fluctuations than low frequency fluc- 
tuations. Emmerling (ref. 8) used an optical method to simul- 
taneously measure pressure fluctuations on a rectangular area. 

His measurements confirmed the previously observed wide variation 
in convection speeds, but he also found that root mean square 
pressure levels were strong functions of the diameter of the 
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pressure sensor. However, Bull and Thomas (ref. 9) suggest that a 
significant portion of the pressure discrepancy is due to the type 
of pressure sensor used in those experiments. Bull and Thomas' 
measurements did confirm a modest increase in rms pressure levels 
with decreasing sensor diameters. Burton (ref. 10) has utilized 
conditional sampling techniques to relate instantaneous wall 
pressure data to the structural features of the turbulent flow. 

He found that "bursts" (see Offen and Kline, ref. 2), which are 
a major source of Reynolds stresses at the wall, are highly 
correlated with the occurrence of large instantaneous, local 
adverse pressure gradients. Burton was unable to link the burst 
structure with the development of the outer structure of the 
turbulent flow . 

The present work is not the first attempt to calculate numeri- 
cally turbulent wall pressure. Deardorff (ref. 11) used a finite 
difference solution to the Navier-Stokes equations for flow in a 
channel to calculate pressure fluctuations. That work has been 
extended by Schumann (ref. 12) for an annular flow. Both were 
able to show good agreement between their computed rms pressure 
levels and experiments. However, because of computational limita- 
tions on spatial resolution and time step, neither was capable of 
simulating energy spectra or space-time correlations. Although 
the Navier-Stokes solution approach is clearly the most funda- 
mentally correct method for calculating turbulent wall pressure 
f luctuations , it is currently incapable of resolving that data 
sufficiently to produce detailed statistical approximations to 
the experimental data. The method employed here attempts to 
simulate the structural features of a turbulent boundary layer 
as observed in experimental measurements, and does not attempt 
to satisfy the equations of motion. Rather, it is anticipated 
that the simulated wall pressure might be useful in conjunction 
with Navier-Stokes solvers by providing a realistic dynamic 
pressure boundary condition. 
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II. SIMULATION STRATEGY 


Ash (ref. 13) has discussed the preliminary strategy employed 
in this work. Essentially, single cycle pressure fluctuations are 
convected over a model surface, and their contribution to the local 
instantaneous pressure field is stored at user- specif ied spatial 
locations with a user-specified time resolution. A Monte Carlo 
method has been used to locate the origin of each pressure fluctu- 
ation event (both spatial location and time), as well as to assign 
frequency and amplitude. The basis for the various distribution 
functions used in the Monte Carlo approach will be discussed here. 

A. Spatial and Temporal Location of 
Individual Pressure Events 

Burton (ref. 10) has shown that the bursts in the near wall 
flow are highly correlated with the wall pressure. Regardless of 
whether the burst structures evolve into the large outer structure 
of the boundary layer or not, the burst data represents a good 
measure of distance and time between individual pressure fluctua- 
tion events. Offen and Kline (ref. 2) indicate that spatial 
separation, Ax, and temporal separation, At, between burst 
events are related by 

Ax = u T At 
or 

u 

Ax = 6* At , (1) 

U oo 

where u is the friction velocity, <$* is the displacement thick- 
ness, is the free stream velocity, and t is dimensionless 

time (t = U^t/6*). 

Burst occurrence measurements from Offen and Kline (ref. 2) 
have been used to construct the probability density function 
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shown in figure 1. From that figure it can be seen that a Gamma 
distribution represents the data reasonably well. Ash (ref. 13) 
has shown that the cumulative probability distribution, constructed 
from the Gamma distribution, is well represented by the Monte 
Carlo distribution function: 

9 JL 

T = 32.2 - - - g 61Q + 72 n 2 + 0.63 tan 2 n (2) 

where n is a uniformly distributed (0 < n < 1) random number. 

Equations (1) and (2) can also be used to assign the distance 
between pressure events, since t is actually Ax — the temporal 
separation between events. 

B. Frequency Associated with Pressure Fluctuation Event 

Bull's (ref. 7) power spectrum measurements have been used as 
a basis for constructing the frequency-generating function. As 
reported by Ash (ref. 13) a cumulative probability function can 
be constructed from the power spectrum. However, the individual, 
single cycle pressure disturbances so generated do not necessarily 
reproduce the prescribed power spectrum. Problems of that type 
encountered in this investigation will be discussed subsequently. 

C. Amplitude of Individual Pressure Events 

Because of the triggering strategy used in Burton’s (ref. 10) 
conditional sampling measurements, it has been possible to infer 
that the distribution of individual pressure event amplitudes is 
Gaussian. Obviously, the mean level is zero and the standard 
deviation is simply the root mean square pressure. 

D. Construction of a Convected Sequence of Pressure Events 

As the individual pressure events are convected downstream, 
their amplitude must decay. Bull (ref. 7) has measured the decay 
rate, and it is well approximated by 
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0 50 100 150 2 

t (U^t/6*) 


Figure 1. Probability distribution of time between bursts. 


(3) 


-4260 v/(u x) 

|p(x + x 0 )|/|p(x 0 )| = 1 - e 

where |p(x + x q )| is the amplitude of the disturbance at a dis- 
tance x from its original location. The above equation does not 
account for frequency- dependent decay. 

Although a substantial amount of evidence indicates that dis- 
turbances are convected downstream with a range of speeds ^0 . 5 £ 

£ 0.85^, a constant and relatively high value of 

u c - 0.8 U m (4) 

has been used here. There appear at this time to be no measurements 
which relate convection speed to size and location of pressure 
disturbances . 

At this point, spatial and temporal spacing, frequency, ampli- 
tude, decay rate, and convection speed models have been discussed. 
The only remaining features are the procedure for constructing an 
orderly progression of pressure events and a determination of the 
required development length for startup of the simulation. 

Obviously, the required model development length can be deter- 
mined using the spatial decay equation. If one assumes that dis- 
turbances which have decayed to an amplitude of one percent of 
their original value are negligible, the required startup, x^, 
is given by 

x D = 0.425 X 10 6 u T /v (5) 

from equation (3). 

The sequence of pressure disturbances is constructed by employ- 
ing equations (1) and (2) along with a random number generator to 
march from the leading edge of the development length to the trail- 
ing edge of the model (region over which pressure is recorded). 

The distance between event birth locations is thus random and 
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distributed in a manner consistent with Often and Kline's (ref. 2) 
burst data. A reference time, t , is frozen during the random 
Ax travel from leading edge to trailing edge, but the time of 
birth of each event, At n > is distributed away from t Q by 


t birth 


t + At 
o n 


( 6 ) 


again, using equation (2) and a random number generator for At . 

Each At n is stored during the spatial travel from leading to 
trailing edge, and an average time step, is computed for 

each complete spatial excursion. A new reference time, t = t Q + 
Atavg 9 is then used > and the procedure is repeated. As each pressure 
event is convected over the desired storage locations, its contri- 
bution is stored at user-specified time intervals. 

Further details of the initial development of this procedure 
are given in Ash, reference 13. At that time, only rms pressure 
level and a pressure signal had been produced. The statistical 
features of the simulation had not been examined, and no attempt 
had been made to optimize the simulation. A detailed study of 
the various factors which influence the quality of the pressure 
simulation has been carried out since that time. 

III. EVALUATION AND OPTIMIZATION OF PRESSURE SIMULATION 


Initially, a single period, sine wave was employed to carry 
the pressure fluctuation event information. The power spectrum 
resulting from the simulation just described is shown in figure 
2. The space-time correlation, defined by 



is shown in figure 3, along with Bull's (ref. 7) measurements. 
Neither of the simulated results was satisfactory. 
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Power spectrum generated by simulation 
from reference 13. 
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Since there was no physical basis for employing a simple sine 
wave as the pressure event carrier, an investigation of other wave 
shapes was undertaken. Four other wave shapes were examined. 

Their names (for reference purposes) and equations are given in 
table 1, along with a sketch of the wave form. A composite power 
spectrum for all five waves (including the original sine wave) is 
shown in figure 4. Except at low frequencies where the spectral 
resolution is not very good (resulting in data scatter), little 
effect of wave shape can be observed on the power spectrum. 

Comparison of the space-time correlations with Bull’s (ref. 7) 
data are shown in figures 5 through 8. It is evident from figure 5 
that the double sine wave affects significantly the space-time 
correlation. Very good agreement in shape and width of positive 
correlation between Bull’s data and the simulation has been achieved. 
At x/6* of 6.28, the correlation is too narrow at its peak, 
but remains in good qualitative agreement. All other wave forms 
(figures 3 and 6 through 8) show varying degrees of agreement. 
Significantly, only the double sine wave carrier caused the space- 
time correlation to possess the negative correlation regions observed 
in Bull’s (ref. 7) measurements. Discrepancies in the correlation 
peaks were not considered significant because they could be adjusted 
by changing the constant in equation (3). 

The wave form investigation was expected to accomplish two 
things. First, it was desired to select a carrier wave which would 
yield space-time correlations which were comparable to the experi- 
mental measurements. Secondly, it was hoped that one or more of 
the wave forms would yield spectral distributions closer to the 
experimental measurements. The second phase was not successful 
per se, but led to two other observations. Not only was the power 
spectrum nearly insensitive to the wave forms examined, but wave 
shapes which were very similar to the double sine wave did not 
produce similar space-time correlations. 

Two important features of the experimental measurements of 
Bull (ref. 7) and Willmarth and Woolridge (ref. 6) had not been 
modeled at this point. Both investigations suggested that 
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Table 1. Carrier wave models 


Name 


Sine 


Double Sine 


Eighth Power 


Jitter 


Skewed 


Equation 


Shape 


Symbol 







Figure 5. Space-time correlation for the double sine 
wave simulation. 
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disturbances were convected with different speeds and that their 
decay was dependent on their characteristic frequencies. As men- 
tioned previously, the constant convection speed may have been 
high, but since no definitive data were available to permit a 
reasonable variable convection speed model, the convection speed 
was not varied. Willmarth and Woolridge (ref. 6) did show the 
influence of disturbance frequency on decay rate, and it was 
decided to investigate the influence of frequency-dependent decay 
on the statistical features of the simulation. 

The filtered space-time correlation decay measurements of 
Willmarth and Woolridge (ref. 6) provide a basis for developing 
a frequency— dependent decay model, but they are not sufficiently 
detailed to permit concise modeling. Their low frequency data was 
taken with a midband Strouhal number of about 0.7 which was more 
than three times higher than the peak frequency. Their high fre- 
quency measurements were centered about a Strouhal number of 5, 
where the energy level was smaller than the peak level by nearly 
two orders of magnitude. In addition, since only two frequency 
bands were employed, there was no way of predicting just how the 
decay rate varied with frequency. 

Decay adjustment with frequency was accomplished by measuring 
the space-time correlation peak in Willmarth and Woolridge* s data 
at x/6* = 2.5. Then, solving the equation 


r pp (2 - 5 ’ 


T ) = 

max 7 


1 - e 


-Av/u x x 


( 8 ) 


for A. Using the two frequency levels it was found that 

Aw ~ 5.6 x 10 6 ^ (9) 

where tu was the midband radian frequency. Because of the form 
of the decay parameter in the computer program, decay was related 
to the disturbance wavelength rather than radian frequency, and 
a reference value for A. The assumed form was 
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|p'(u>, x)| 

= 1 - e 

I p ’ ( 0) , o)| 


( 10 ) 


-3230 vX/(xu T X r ) 


where A r is the wavelength of a reference disturbance, given by 

A r = 0.8 6* y (11) 

with y a parameter (taken as 30 in the test) used to relate the 
reference disturbance length to a small disturbance length (0.8 6*). 

A simulation was generated using the linearly dependent fre- 
quency (wavelength) decay function given in equation (10). The 
power spectrum resulting from that simulation is shown in figure 9. 

It can be seen that frequency— dependent decay has significantly 
altered the power spectrum. Furthermore, the alteration — particularly 
at high frequencies — is favorable. Unfortunately, examination of 
the associated space-time correlation, shown in figure 10, indicates 
that linearly varying decay has had a devastating effect. 

Three significant observations from this study can be made at 
this point: (1) the double sine wave carrier is clearly superior 

to the other wave forms examined with regard to producing reason- 
able space-time correlation curves; (2) wave form does not affect 
significantly the power spectrum, and the Monte Carlo frequency 
generating function used thus far does not appear to be acceptable 
in simulating the power spectrum because of the absence of a 
clearly defined peak in the energy spectrum; and (3) frequency- 
dependent decay is a major influence on the shape of both the 
power spectrum and the space-time correlation curves. These 
observations permitted a threefold adjustment of the simulation 
program. Each adjustment will be addressed separately. 

A. Wave Form 


The double sine wave disturbance carrier was used for the 
remainder of the simulations. Its obvious higher harmonic contri- 
bution to the wall pressure energy spectrum did not signif icantly 
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Frequency (Hz) 


Figure 9. Wall pressure power spectrum for sinusoidal disturb- 
ances, decay varying linearly with frequency. 
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Figure 10. Space-time correlation for sinusoidal waves with 
decay a linear function of frequency. 


alter that spectrum'. The fact that different waves nominally with 
the same shape were not capable of preserving the desired features 
of the space-time correlation suggests that the higher harmonic 
aspect of that wave carrier may have physical significance. 
Regardless of the physical significance, no further wave form 
studies appeared warranted. 


B. Frequency-Generating Function 

An adjustment was certainly required for the Monte Carlo fre- 
quency generator. However, it was uncertain at this point just 
how the combined effects of frequency-dependent decay and a modi- 
fied frequency generator might interact. Rather than examine the 
effects separately, the composite spectrum shown in figure 4 was 
selected as a guide for adjusting the frequency generator and the 
decay variation, discussed subsequently, was modified simultaneously. 

A curve was faired through the simulated spectral data shown 
in figure 4 and a ratio of desired level (based on Bull's data) 
to generated level was established. The resulting spectrum, 
gotten by multiplying the ratio of Bull's values to the simulated 
values by Bull's original spectrum, is shown in figure 11. Using 
the resulting spectrum, the cumulative probability function was 
obtained by graphically integrating the area under the power 
spectrum curve. That distribution is also shown in figure 11. 

Using the features of the cumulative probability distribution 
function (particularly the upper and lower frequency limit asymp- 
totes) shown in figure 11, the dimensionless frequency-generating 
function given by 


a) = 0.523 


- i 


+ 0.799 n 2 / 3 - 0.785 n 


(1 - n) 


74 


( 12 ) 


was constructed. No significant deviation between the curve of 
equation (12) and the curve shown in figure 11 was observed. 
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Figure 11. Modified power spectrum and cumulative probability 
required for frequency generator. 


24 


Cumulatjive Probability 



C. Frequency-Dependent Decay 

Based on the drastic effect of frequency- dependent decay on 
the space-time correlation curves, as shown in figure 10, the 
linear variation of decay with frequency was abandoned. The 
difficulty was due to the fact that at the lowest frequencies 
decay rates were two orders of magnitude slower than at the high- 
est frequencies — causing an inappropriate overemphasis on the low 
frequency waves and resulting in the greatly widened positive 
space-time correlation regions shown in figure 10. 

Frequency dependent decay was still desirable from the experi- 
mental implications, as well as its favorable effect on the power 
spectrum. An exponentially varying frequency dependent decay 
function was selected because it was capable of maintaining nearly 
constant decay rates at the lower frequencies (below Strouhal 
numbers of about 0.2), while increasing significantly the rate of 
decay at the higher frequencies. The form used was given by 

|p(w, x)| -3000 v f((jj)/(u x) 

= 1 - e T (13 

|p(o), 0)| 


where f(a)) was given by 


(to/a3 o ) -(X r /X) 

(gd) = e 


(14) 


with a) Q specified in the same manner as X p in the linear case 
(y was taken as 10 in this case). 

The power spectrum and space-time correlations generated by a 
simulation using the modifications just described are shown in 
figures 12 and 13 respectively. Good agreement between the simu- 
lation and experimental data is evident. The deviation in decay 
rate indicated in figure 13 is small (about 10 percent) and can 
probably be reduced further by adjustments in the decay function. 
However, the agreement was considered adequate at this point. 
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Figure 12. Power spectrum after optimization. 





In concluding this section, it should be noted that no formal 
discussion of rms pressure level has been made thus far in this 
report. It is extremely simple to adjust rms pressure level in 
this simulation because one need only adjust the prefix of the 
Gaussian amplitude generator. The amplitude of the rms pressure 
level has no affect on the other statistical properties of the 
simulation. The desired rms pressure level was assumed to be 
3 t w — slightly higher than either Bull’s (ref. 7) or Willmarth and 
Woolridge’s (ref. 6) measurements, but consistent with the observa- 
tions of Emmerling (ref. 8) and Bull and Thomas (ref. 9). As 
evidenced by the comment statement in the computer program listed 
in the appendix of the present report, a somewhat reduced prefix 
magnitude (20 percent reduction) was required to yield a nominal 

p of 3 t at the five simulation stations. Variations of 

*rms w 

each of the rms pressure levels from the specified value were 
about plus or minus five percent. 

IV. RESULTS AND DISCUSSION 

This investigation has shown that it is possible to simulate 
turbulent wall pressure fluctuations using a Monte Carlo approach. 
The simulation is capable of satisfactorily reproducing rms pressure 
levels, power spectra, and space-time correlations (in the direc- 
tion of flow). The simulation permits arbitrary resolution of 
pressure data in both space and time, making the output directly 
accessible to other computational analyses. 

The latest version of the FORTRAN computer code for generating 
the simulation is included in the appendix to this report. The 
space-time correlation calculations generated at the end of the 
listed simulation program are not correct, but were merely used 
to check the computation prior to storage on tape. The simulation 
of 0.767 seconds of pressure data at five spatial locations with a 
resolution of 46.5 microseconds required nearly 2800 seconds of 
CDC 6600 computer time. 

Obviously, a great deal of computer time is required to gener- 
ate a simulation. However, generation of a 1 or 2 second simulation 
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can be used to produce an essentially infinite time simulation simply 
by randomly accessing segments of the simulation and piecing them 
together. Since the data is on tape, the simulation output can 
be read, one time step at a time, into the other computer programs, 
thereby alleviating the storage requirements for further analyses. 

Arbitrary spatial resolution is also limited by the storage 
capacity of the computer, i.e., pressure data cannot be stored at 
an infinite number of spatial locations. However, that problem 
can also be bypassed by taking advantage of the fact that the 
random number generator in a computer reproduces itself. That is, 
if the same startup number (seed) is used for all runs and the 
length of the model does not change (the last storage location 
must be fixed), the simulation can be repeated. Since the space 
and time resolution of the output is arbitrary, the computer will 
go through an identical chain of random numbers and, by changing 
the storage locations for the output, an arbitrary number of 
pressure histories for the specified spatial locations can be 
generated and stored on one or more tapes. 

Besides the ability to generate infinite duration one-dimensional 
pressure histories, a two-dimensional simulation can be produced in 
much the same manner. If the transverse spacing, Ay, between 
storage locations is greater than about 25 6*, the same simulation 
can be randomly accessed to produce pressure histories for succes- 
sive strips, since R (x, Ay, t) would be nominally zero. If the 
spacing is less than 25 6*, a similar approach can be used. 

When two-dimensional simulations are required where Ay < 25 6*. 
there are a countable number of locations that occupy a strip 25 6* 
wide. Using the procedure above to lay down histories on the loca- 
tion lines that exceed 25 6*, the intermediate strips can be gener- 
ated by alternately assigning data from the independent strips. 

The transverse space-time correlation could be used to construct 
a generating function to yield such a simulation, but further 
research would be required to augment that case. 

Finally, since spatial data scales with 6*, and time data 
with 5*/U ot , a given simulation can be used for more than one 


29 



flow condition, but the user loses control of space and time reso- 
lution of the output, as well as the skin friction coefficient 
(since u^VU^ has already been set in the simulation). The 
usefulness of a general simulation for all turbulent boundary 
layers is therefore rather restricted. 


The author would like to thank Ramakrishna Balasubramanian 
for his assistance in constructing the computer codes used in 
this investigation. 



APPENDIX 


PRESSURE SIMULATION FORTRAN PROGRAM 


on n n nnnnno 


JOB « 1 • 7700 • 3000 00 * 1 0000 • A4677 R4313 101409 BIN34 

USLR.UALA.RAMAKR 1SHNAN 000029300E 37200 NAS OPU 

L INECNT (50000 ) 

REQUEST . T APE 1 * HY • SA V TP . R I L » RSB • PRF 1 

RE W I NO ( T APE I ) 

RUNCS ) 

LGO. 

PROGRAM S I MU ( INPUT • OUT PUT . TAPE5= INPUT * TAPE6 = OUTPUT * T APE 1 • TAPF2 ) 
DIMENSION P ( 5 * 1 6500 ) 

DIMENSION SUMX (8000 ) •SUMY ( 0000 ) 

DIMENSION PA(5).PR(5) 

EQUIVALENCE (P (66000 ) .SUMX( 1 ) ) * (P ( 74000) . SUM / < 1 ) ) 

NR=5SNW=6 

C ND I M = M A X I MUM TIME DIMENSION 

NDIM= 1 6500 
CD I M-ND1 M 

CORX IS THE ADJUSTMENT FOR DX STREAK. 

IF MORE THAN ONE BURST PER STREAK • CORX IS NOT UNITY 
CORX= l • 

CORT IS THE TIME ADJUSTMENT FOR DT-STRFAK 
CORT = 1 . 

XD=DEVELOPMENT LENGTH (M) 

XM=MODEL LENGTH (M) 

PAR IS THE DECAY ADJUSTMENT ACCOUNTING FOR DISTURBANCE S I ZF 
PAR= 1 0 » 

XD=4 • 1 45 
I MJ = 0 
XM=0 • 0 
DXT=0. 008 
157 XD=XD+XM 
XT = XD+XM 

DXT = NODE SPACING ON MODEL (M) 

NXM = 5 

XSMX=XD+DXT#NXM 

XSMX IS THE RIGHT-MOST LOCATION STORING PRESSURE 
US=33.5iBLT=0.0408$UFRIC= 1 • 2587 SCNU= 0 • 0000 1 395 
RHO = 1 .2 
DPT=BLT/B 

ENTER OTHER DISPLACEMENT THICKNESS HERE IF DESIRED (CM) 
CALCULATION OF WALL SHEAR STRESS. TW 
TW=RH0*UFRIC*UFRI C 
C CALCULATION OF RMS PRESSURE FLUCTUATION 

PRMS-3 • * T W 



o n 


CALCULATION OF MINIMUM RELEVANT DISTURBANCE LENGTH 
TWV IS TIME DURATION OF DISTURBANCE 
TWV=DPT/US 
C LENGTH OF SMALL DISTURBANCE 

CLM=0.8*US*TWV 
C DECAY BIAS PARAMETER 

CWV=PAR»CLM 

C MISCELLANEOUS CONSTANTS 

P I “3. 1 4 1 5926 
TPI =2.*P1 
HP I =P 1/2. 

C0=?.5l 551 7 
Cl =0 • 802853 
C2 = 0 • 0 1 0320 
D 1 =1 • 432700 
02 = 0 • 109269 
D3 = 0 • 00 1 300 

C STARTER FOR RANDOM NUMBER GENERATOR 

XSTART = 77653 • 

RNM=UR AN (XSTART ) 

XST ART = 0*0 

C CALCULATION OF NOMINAL PEAK FREQUENCY 

FPEAK = 0 *20574 *US/ ( TP I *DPT ) 

FMAX=1 0 .*FPEAK 
C COMPATABLE TIME STEP 

DTT = | « /FMAX/ 1 0 • 

C A USER SPECIFIED TIME STEP CAN BE ENTERED HERE 

C THIS TIME STEP lS THE TIME INCREMENT USED IN THE RESOLUTION OF THE 

C OUTPUT THE INTERNAL FLUCTUATION TIME STEP IS RANDOM 

C MAXIMUM TIME IS CONSTRAINED BY COMPUTER STORAGE • SET THE NUMBER 

C OF ALLOWABLE TIMESTEPS NTM 

NTM=16500 

C INITIAL TIME IS TO (SEC) 

TO = 0. 

DT A y/G“ O • 

C PREFIXES FOR RANDOM NUMBER CALCULATIONS 

PX = DPT*UFR I C/US 
PT=DPT/US 
PW=US/DPT 

C CALCULATION OF REQUIRED START UP TIME FOR SIMULATION 

TSO=l • 0*XT /US 

nmin=tso/dtt 

TMAX1 =NMIN*DTT 

C IF NM IN IS GREATER THAN NTM* NTM IS OVERRIDDEN 



IF<NMIN-NTM)201 *201 *202 
202 NTM=NM 1 N 
201 CONTINUE 

TMAX=NT M*DTT 
WRITE(6*I00) DTT,TMAX 

1 00 FORMAT <5X* *DT = * *F1 O .0* 1 Dx . *TMAX = -» »F 1 O.A , // ) 

TMAX5=CD IM*0TT 
TMAX=TMAX+TMAX I 
NTM2=NDI M'NM I N + 2 
T REF "TMAX 1 
NREL=0 

C INITIALIZE PRESSURE ARRAY 

DO 52 NX=1 *NXM 
DO 52 NT= 1 « ND I M 
52 P ( NX • NT )=0.0 

nflg=o 

TSUB=0. 

C INITIALIZE LOCATION AND TIME BASE* ETC. 

1 X=0. 

NCT = 0 
DTSUM=0 . 

T 0= TO + DTAVG 
IF (T0-TMAX1 >2,150,150 

150 IF ( NFL G > 151*151*15 

151 NFL G= 1 
TMAX] =TMAX 

DO 154 1=1 ,NXM 

DO 152 J=NM I N • ND I M 
JJ=J+1-NMIN 
P( I ,JJ)=P( I ,J) 

152 CONTINUE 

DO 1 53 J= N TM2 , ND I M 
P ( I . J ) =0.0 

153 CONTINUE 

154 CONTINUE 
TSUE3 = TREF 

T MAXS = T MA XS + TREF 
NREL=NM I N 

2 RNM=URAN IXSTART > 

RNM=0. 005+0. 99*RNM 

C CALCULATION OF DX USING RANDOM NUMBER RNM 

HP I I = HP I *RNM 

DX=PX*(32.2-2/(RNM + 0.619) + 72.*RNM**2 + 0.63*TAN(HPI I > ) 
OX=CORX*DX 
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X - X+DX 

«NM*URAN ( XSTART ) 

RNM=0 • 005+0. 99*RNM 

C CALCULATION OF RADIAN FREQ* FROM NEW RNM 

SRNM*RNM**0 *6667 
RRNM* C 1 *-RNM )**0.74 
FRNM» I /RRNM- I 

W=PW* (0*523*FRNM+0*799*SRNM-0«7e5*RNM ) 

F*W/TPI 

TP=l/F 

DXE=0.8*U5*TP 
XO r X—DXE 

C XO IS THE ORIGIN OF THE SINE WAVE FLUCTUATION 

C X IS THE FRONT OF THE SINE WAVE 

C XS IS THE FIRST STATION AT WHICH P IS RECORDED* 

C XS AND NX I WILL BE TAKEN AS THE FIRST oT AT I ON VALUES THEN OVER I DDEN 

XS=XD 
NX I = 1 

C CHECK TO SEE IF THE DISTURBANCE IS OVER THE MODEL 

IF(X-XD)5*5.3 

IF THE DISTURBANCE IS OVER THE MODEL * HAS IT PASSED THE LAST DATA 
STAT I ON 

3 IF CX-XSMX )4„t , 1 

C NX I IS THE NUMBER OF THE FIRST STORAGE LOCATION 

4 NX I = (X-XD J/DXT+l .9999 
CXNI =NXI 

XS r XD+CXN I *DXT 

5 CONTINUE 

C GENERATION OF RANDOM TIME STEP 

6 RNM=uRAN ( XSTART ) 

RNM = 0* 005+0 • 99*RNM 
HP I T = HP I *RNM 

DT*PT* C 32.2-2*/ (RNM+O. 61 9 ) + 72 • *RNM**2 + 0 • 63* T AN ( HP I I ) ) 

DT=CORT*DT 
T * TO+DT 
NCT=NCT+1 
DTSUM=DTSUM+DT 

dtavg=dtsum/nct 

C GENERATION OF GAUSSIAN RANDOM PRESSURE AMPLITUDE 

RNM=URAN (XSTART ) 

C IND=RNM+0*5 
IND=C1ND 
C I ND= I ND 

PPP 3 2 * * ( 1 * — C I ND ) — 1 • 
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ARGR=RNM/(1 ,+ClNO) 

arg=i •/ < argr* argr ) 

C T = ALOG ( ARG ) 

CM=SORT <CT ) 

PMG = CM- (CO+CM* ( Cl +CM*C2 ) )/ ( H CM* <D1 + CM* (D2+CM*D3 > ) ) 

XS=XS-OXT 

PE;=PRMS*PMG*PPP 

C MODIFICATION OF PRESSURE AMPLITUDE" lO ADJUST RMS PRESSURE LEVEL 

PE=0 • 8*PE 

C DO LOOP FOR STEPPING THROUGH MODEL STORAGE LOCATIONS 

DO 1 A NX = NX I « NX M 
C MODEL STATION X-LOCATION 

XS = XS+DXT 

C ARRIVAL TIME OF PRESSURE FLUCTUATION 

DXS= XS— X 

T GO= 1 .25*DXS/US + T 
C FLUCTUATION DEPARTURE TIME 

D XO= XS-XO 
TSTP- 1 . 25*DX0/US+T 
C DOES TSTP EXCEED T MAX 

IF (TSTP-TMAXS ) 9 « 9 , 7 
7 IF (TGO-TMAXS)B.B, 13 
B NSTOP=NDIM 
GO TO 10 
9 NSTOP=TSTP/DTT 
JO NGO=TGO/DTT 

I F (NGO-NREL ) 99,99*98 
9B CONTINUE 

C DO LOOP FOR SUCCESSIVE TIME CONTRIBUTIONS TO THE SAME X LOCATION 

DO 12 NT=NGO«NSTOP 
T C = NT *DT T 

THET=TP I * (TC-TGO }/TP 

DELT=TC-T 

XO T = 0 • 8 ♦US* DEL T 

IF (XOT-0.0OO5 119,19,18 

18 CONTINUE 

ARGX=-3000*CNU/ (XOT*UFR I C ) 

C FREQUENCY DEPENDENT VARIABLE DECAY RA I t AQ JUS I MtN I 

ARGAR=— CWV/DXE 
ARGX=ARGX*EXP ( ARGAR ) 

DECA= 1 —EXP ( A RGX ) 

GO TO 23 

19 CONTINUE 
DECA = 1 • 



23 CONTINUE 

DP=PE*5 IN ( THET ) * ( ] • -COS ( THET ) ) 

DP = DP* DEC A 
NT IME = NT-NREL 

P ( NX * NT I ME )=P (NX, NT I ME > +DP 
12 CONTINUE 
99 CONTINUE 

13 CONTINUE 

14 CONTINUE 

130 FORMAT (5X,F10.6*6E1 4 ,6 ) 

GO TO 2 
15 CONTINUE 

NWRT=NTM/8 
CNT = NTM 
DO 11 J=li5 
PA ( J ) =0*0 
11 PR ( J ) = 0 • O 

DO 21 J= 1 «5 
DO 20 1=1 * NT M 

PA(J)=Pfl(J)+P(J.| )/CNT 
PR ( J ) =PR ( J ) + IP ( J« 1 )**2/CNT ) 

20 CONTINUE 

PR( J)=PR( J)**0.5 

21 CONTINUE 

DO 22 J=l «5 
DO 22 I = I • NTM 

22 PI Jt I )=PC Jt I )-PA< J) 

NYM= 1 OOOSNZM = 200 1 
ANN=200 1 

DO 24 KK= 1 , 4 
DO 25 K=1 ,NYM 
SUM=0,tSUMl =0. 

K I =K— 1 
NZ=NZM-Kl 
DO 40 M=1 ,NZ 

SUM=SUM+P ( 1 • M ) *P ( 1 ♦ M+K | ) /PR (1 >**2 
SUM 1 =SUM I +P ( 1 »M)*P(KK+I • M + K 1 )/P« ( 1 ) /PR (KK + J ) 
40 CONTINUE 

SUMXIK 1 =SUM/ANN 
SUMY (K ) =SUM1 /ANN 
25 CONTINUE 

1 F (KK.GT. 1 ) GO TO 75 
WR ITE CNWt 191 1 
191 FORMAT C 1 HI ) 



WRITE<NW*200 ) <SUMX<K ) « K= 1 * 1000) 

200 FORMAT (5E20 .6 ) 

75 CONTINUE. 

WRITF (NW* 1 91 ) 

WR I TE < NW.200 ) ( SUMY < K )iK= I « 1 000 ) 

NXX = 1 000 
DO 80 K= 1 • NXX 
SUM=0.»SUM| =0* 

K 1 =K-1 

DO 90 M=K • NYM 

SUM = SUM + P ( 1 *M >*P ( 1 ,M-K 1 )/PW ( 1 > **2 
SUM 1 = SUM 1 +P ( 1 *M )*P{KK+1 i M— K 1 )/PR ( 1 )/PR(KK+l 
90 CONTINUE 

SUMXCK ) = SUM/ ANN 
SUMYtK )=SUM1/ANN 
80 CONTINUE 

I F ( KK • GT • 1 ) GO TO 793 
WRITEINW, 191 ) 

WRITE (NW* 200 ) (SUMX(K ) *K= 1 . 1 000 ) 

793 CONTINUE 

WR I TE ( NW • 191 ) 

WRITE (NW.200 ) (SUMY(K > * K ~ 1 , 1 OOO ) 

24 CONTINUE 
JXY=5500 
DO 93 I JJ=1 *3 
J1 =( I JJ-1 ) * JXY + 1 
J2= I JJ* JXY 
DO 94 KK=1 *5 

WR ITE ( 1 I(P(KK»JL)tJL=JI «J2) 

ENDFILE 1 
94 CONTINUE 
93 CONTINUE 
STOP 
END 
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